#  This is part of the Scopia data visualization example collection
#
#
#  Copyright (C) 2008  Scopia Visual Interfaces Systems 
#                      http://www.scopia.es
#
#  author: Alejandro Conty
#  contributors:
#
#  This program is free software: you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation, either version 3 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
from numpy import *

def buildMultipleBezierCurves( pointMatrix, steps = 8 ):
    """Build an array of bezier curves (sequences of points) from an array
       of sets of control points that define each one. pointMatrix is a matrix 
       of the form:
        [ [p1, p2, p3, p4],
          [p1, p2, p3, p4],
          ... ]
       With a shape of (N, 4, 3).
       And the return array will have shape (N, steps+1, 3)"""
    # t variable and one minus t variable
    t = arange( 0, 1.0 + 0.5/steps, 1.0/steps, float32)
    omt = 1.0 - t
    # Bezier matrix for the given number of steps
    M = zeros( (steps+1, 4), float32 )
    M[:,0] = omt*omt*omt
    M[:,1] = 3*omt*omt*t
    M[:,2] = 3*omt*t*t
    M[:,3] = t*t*t
    # product of the matrix by the control points of every curve in the input array
    curves = sum(M[..., newaxis] * pointMatrix[:, newaxis, ...], axis=2)
    return curves

def get2DNormalsFromCurves( curves ):
    """Given an array of curves (sequences of points), calculate some kind
       of "normal" vectors, perpendicular to the curve at every point in
       the sequence (laying in XY plane). We need this for building quads."""
    mdeltas = zeros( curves.shape, float32 )
    # first all the vectors out of every point minus the previous one
    mdeltas[:,:-1] = curves[:,1:] - curves[:,:-1]
    # last point has the same normal as the previous one
    mdeltas[:,-1] = mdeltas[:,-2]
    # then swap X and Y and change one sign
    temp = mdeltas[:,:,1].copy()
    mdeltas[:,:,1] = - mdeltas[:,:,0]
    mdeltas[:,:,0] = temp
    # normalize all of them
    lengths = sqrt(sum( mdeltas*mdeltas, axis=2 ))
    mdeltas /= lengths[...,newaxis]
    return mdeltas

def getQuadsFromMultipleCurves( curves, normals, widths ):
    """Having several curves and the corresponding "normals" for every point
       build quad lists following the trayectories with the given width"""
    quads = zeros( (curves.shape[0], curves.shape[1] - 1, 4, 3), float32 )
    # unoptimized code :
    #quads[:,:,0] = curves[:,:-1] - normals[:,:-1] * widths[:,newaxis,newaxis]
    #quads[:,:,1] = curves[:,:-1] + normals[:,:-1] * widths[:,newaxis,newaxis]
    #quads[:,:,2] = curves[:,1:] + normals[:,1:] * widths[:,newaxis,newaxis]
    #quads[:,:,3] = curves[:,1:] - normals[:,1:] * widths[:,newaxis,newaxis]
    # optimized code :
    widths = widths[:,newaxis,newaxis]
    quads[:,:,3] = quads[:,:,0] = - normals[:,:-1]
    quads[:,:,1] = quads[:,:,2] = normals[:,1:]
    quads[:,:,0] *= widths
    quads[:,:,0] += curves[:,:-1]
    quads[:,:,1] *= widths
    quads[:,:,1] += curves[:,:-1]
    quads[:,:,2] *= widths
    quads[:,:,2] += curves[:,1:]
    quads[:,:,3] *= widths
    quads[:,:,3] += curves[:,1:]
    # we don't need the quads to be grouped in 4 points sets, just a sequence
    # of points as opengl expect them. That's why we "eat" one dimension here
    quads.shape = (curves.shape[0], (curves.shape[1] - 1) * 4, 3)
    return quads

def getMultipleBezierQuads( pointMatrix, width, steps = 8 ):
    """This is the useful function in this module. It already
       returns bezier curves made out of quads from the array
       pointMatrix. Which is a sequence of sets of control points
       defining a bezier. Width is used for building the quads.
       pointMatrix is a matrix of the form:
        [ [p1, p2, p3, p4],
          [p1, p2, p3, p4],
          ... ]
       With a shape of (N, 4, 3).
       And the return array will have shape (N, steps*4, 3)"""
    curves = buildMultipleBezierCurves( pointMatrix, steps )
    normals = get2DNormalsFromCurves( curves )
    quads = getQuadsFromMultipleCurves( curves, normals, width )
    return quads

